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By shifting the reference system for the local-density approximation (LDA) from the electron gas 
to other model systems one obtains a new class of density functionals, which by design account for 
the correlations present in the chosen reference system. This strategy is illustrated by constructing 
an explicit LDA for the one-dimensional Hubbard model. While the traditional ab initio LDA is 
based on a Fermi liquid (the electron gas), this one is based on a Luttinger liquid. First applications 
to inhomogeneous Hubbard models, including one containing a localized impurity, are reported. 
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Density-functional theory (DFT) [1] is the basis of al- 
most all of todays electronic-structure theory, and much 
of materials science and quantum chemistry. Many-body 
effects enter DFT via the exchange-correlation (xc) func- 
tional, which is commonly approximated by the local- 
density approximation (LDA) [1]. The essence of the 
LDA is to locally approximate the xc energy of the in- 
homogeneous system under study by that of the homo- 
geneous electron gas. This electron gas plays the role 
of a reference system, whose correlations are transfered 
by the LDA into the DFT description of the inhomoge- 
neous system. The most popular improvement upon the 
LDA are generalized gradient approximations [2] , whose 
basic philosophy is to abandon the requirement of homo- 
geneity of the reference system. This system, however, is 
normally still the interacting electron gas [2] . 

In the present paper we propose to explore a different 
paradigm for the construction of novel density function- 
als: instead of sticking to the electron g reference 
system, and abandoning homogeneity, it may sometimes 
be advantageous to do the reverse: stick to homogeneity 
(and thus to the LDA) but abandon the electron gas as 
a reference system. The new reference system is chosen 
such that it accounts for the correlations present in the 
inhomogeneous system under study. 

The only requirement for the reference system is that 
in the absence of any spatially varying external potential 
its xc energy must be known exactly or to a high degree of 
numerical precision. Besides the electron gas (or Jellium 
model) there are many other physically interesting model 
systems that satisfy this criterium. Most notably among 
these is a large class of low-dimensional models which 
can be solved exactly by Bethe Ansatz (BA) techniques 
or bosonisation (in one dimension, e.g., the repulsive and 
the attractive Hubbard model, the hard-core Fermi and 



Bose gases, the Heisenberg, the supersymmetric t-J, and 
the Tomonaga-Luttinger model [3,4]). The solutions to 
these models in the homogeneous case can be used in- 
stead of the electron gas to construct LDA functionals 
that can then be applied to study these models also in in- 
homogeneous situations. The main advantage offered by 
a DFT treatment of such models is the gain in simplicity 
that arises from mapping the inhomogeneous interacting 
many-body system onto a nonintcracting auxiliary sys- 
tem, which is diagonalized much more easily. 

Below we implement these ideas for the one- 
dimensional Hubbard model (1DHM), by constructing an 
LDA based on the exact Bethe Ansatz solution of Lieb 
and Wu [5]. A DFT treatment of the Hubbard model 
has been pioneered by Gunnarsson and Schonhammer in 
Ref. [6], but the LDA-type functional they proposed has 
in practice often led to disapointing results [7], and was 
criticized as not being a proper LDA since it was not 
based on the exact solution of a homogeneous reference 
system [8]. As a consequence, more complicated approx- 
imation schemes, such as self-interaction corrections, are 
often employed [9]. An attempt to base a proper LDA 
for the Hubbard model on the BA was made in Ref. [10], 
but the formulation of that work has not been widely 
applied, probably because no explicit expression for the 
resulting xc functional was provided. 

Motivated by our above analysis of the possibility of 
switching reference systems for the LDA we construct, 
in the present paper, an explicit and simple BA xc func- 
tional and apply it to a variety of inhomogeneous Hub- 
bard models, among them an impurity model for previ- 
ously unattainable system sizes. This functional, denoted 
the BA-LDA, has built into it the Luttinger-liquid cor- 
relations present in the 1DHM [3,4], in the same way in 
which the conventional LDA has built into it the Fermi- 
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liquid correlations present in the electron gas. 

The Hamiltonian for the homogeneous 1DHM is, in 
standard notation, 



iT c iT c a c a . 



(1) 



Here t is the kinetic energy (in the following taken to 
be the unit of energy) and U the interaction (considered 
a fixed parameter, characterizing the Hamiltonian). To 
construct an LDA we first develop a parametrization for 
the total energy per site, as a function of U and n (the 
filling factor, a constant in the homogeneous case). Our 
parametrization interpolates analytically between three 
limiting cases in which explicit results can be extracted 
from the Bethe-Ansatz solution [3,5]: (i) U — ► and any 
n < 1, (ii) U — > oo and any n < 1, and (iii) n = 1 and 
any U [11]. The expression 
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where /3(f) is an (n independent) number which is de- 
termined for any given value of U from 
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and Jo and J\ are zero and first order Bessel functions, 
recovers all three limits: The right-hand side of Eq. (3) 
is the exact BA expression for the total energy at half 
filling. The parameter (3 is thus determined such that at 
n = 1 Eq. (2) becomes exact. On the other hand, Eq. (2) 
is already of the algebraic form of the exact results for 
the limits U = and U — > oo, in which (3 — 2 and [3=1, 
respectively. In these limits the integral in Eq. (3) can 
be calculated analytically, and one indeed recovers these 
values for (3{U) [12]. 

Eq. (2) with (3(U) determined from (3) is thus exact in 
the three situations mentioned above. In order to check 
whether these equations provide a reasonable approxi- 
mation also between these limits, we have numerically 
solved the Lieb-Wu integral equations for the full Bethe 
Ansatz solution [5] and compared the resulting total and 
xc energies with the one obtained from Eqs. (2) and (3). 
We find that both agree to within at most a few percent, 
even for values of the parameters far away from the exact 
limiting cases. This is illustrated in the inset of Fig. 1. 

In order to extract the exchange-correlation energy 
ef c A (n, U) from the total-energy expression (2), one fol- 
lows the usual prescription of DFT [1] and subtracts the 
Hartree energy and the noninteracting kinetic energy. 
This latter energy is simply given by substituting /3 = 2 
in Eq. (2) . In the spirit of the usual electron gas LDA we 
can then construct an LDA for the 1DHM, i.e., 
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BA-LDA 
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where n, = '^2 a {c\ (T Ci C! ) . Given this expression for the xc 
functional, ground-state properties of Hubbard models 
subject to a wide spectrum of inhomogeneities can be 
calculated from DFT. Note that in such a calculation 
Eq. (3) must be solved only once for any given value of 
U, i.e., determination of (3 takes place outside the self- 
consistency cycle of DFT. 

As a first numerical example we apply the BA-LDA to 
a finite 1DHM with open boundary conditions, and cal- 
culate its ground-state energy as a function of the number 
of sites L. In Fig. 1 we compare our results with those ob- 
tained from exact (Lanczos) diagonalization of the same 
system. We see that around L = 7 there is a crossover 
between exact and approximate data points (separately 
for even and odd values of L). After that, the deviation 
between both sets of data saturates to a value near the 
intrinsic error of the interpolation (2), indicated by the 
error bar. As an example for a truly inhomogeneous sys- 
tem we now add an on-site potential ^2 i(7 Vic\ a Ci a to the 
Hamiltonian. Our results for a binary potential (with 
Vi = — 1 on the odd sites and Vi = +1 on the even ones) 
are displayed and compared with exact diagonalization 
in Fig. 2. 

We have performed similar calculations also for sev- 
eral other values of U and other external potentials m. 
Our conclusions from these calculations are: (i) The ac- 
curacy of the BA-LDA total energy is typically of the 
order of a few percent, and much better than that near 
crossovers and near the limits at which equation (2) ex- 
actly represents the underlying homogeneous reference 
data (U = 0, U — > oo, n = 1). (ii) Unlike traditional 
methods, the quality of the BA-LDA does not deteriorate 
as L increases, and the computational effort associated 
with it is that of diagonalizing a noninteracting system. 
Fully self-consistent calculations for systems with tens of 
thousands of sites are thus possible. This is a unique, and 
rather desirable, feature of the BA-DFT, as compared to 
traditional methods. 

Interestingly, when one of the site occupation numbers 
comes very close to 1 (typically within less than 5 x 10 -3 ), 
the self-consistency cycle associated with the Hubbard 
model Kohn-Sham equations docs not necessarily con- 
verge. In the homogeneous case n = 1 (half filling) marks 
the Mott metal-insulator transition associated with the 
opening of a gap in the energy spectrum [3-5] . Whenever 
in a metallic system one of the site occupation numbers 
comes close to 1, a local approximation, such as the BA- 
LDA, thus treats the system at that site as if it were an 
insulator, in spite of the fact that the metallic (Luttinger 
liquid) correlations are very different from those of the 
Mott insulator. The 1DHM thus constitutes a theoret- 
ical laboratory in which the band-gap problem of DFT 
can be studied [6,10]. Results obtained with the present 
functional will be reported in a forthcoming publication. 

After these preliminary investigations we now consider 
a case that illustrates the full power of the BA-LDA ap- 
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proach: a localized impurity. For this kind of system tra- 
ditional approaches face the problem of slow convergence 
to the thermodynamic limit. Different types of impurities 
in the Hubbard model have been studied in the literature 
by various techniques [6,13]. Here we model the impu- 
rity by choosing vi = — 1 at the impurity site and v = 
everywhere else, so that electrons will be dragged to the 
impurity site. The density distribution for the 201-site 
system is displayed in Fig. 3, while the convergence to 
the thermodynamic limit is illustrated in Fig. 4, which 
also contains results obtained for a much more attractive 
impurity with vj = —10. 

The difference between open and periodic boundary 
conditions becomes small only when the system size L 
exceeds the damping length of the Friedel oscillations 
originating at the surface in the open case. To bring 
out clearly the effect of the impurity, regardless of the 
choice of boundary conditions, we have, in both figures, 
subtracted the results for the same boundary conditions 
without the impurity. What remains are the Friedel os- 
cillations originating at the impurity. These oscillations 
on their own significantly slow down the convergence 
of the total energy to the thermodynamic limit: Fig. 4 
shows that, as expected in the thermodynamic limit, the 
impurity-added contribution to the total energy scales 
linearly with the impurity concentration 1/L. The im- 
purity problem thus illustrates an area in which the BA- 
LDA can be useful, since systems of the size required to 
approximate the thermodynamic limit to within a per- 
cent or better are hard to study with traditional meth- 
ods, in particular for periodic boundary conditions. 

For small systems, where exact diagonalization is pos- 
sible, we also compared the density distributions ob- 
tained from the BA-LDA with the exact ones. Both agree 
quantitatively. For larger systems one can use the BA- 
LDA results to study the asymptotic algebraic decay of 
the oscillations. For the case depicted in Fig. 3 we find, 
for example, that the oscillations decay as l/x 7 , where x 
is the distance to the impurity site and j(U = 6) = 1.30. 
This exponent is a nonuniversal (interaction-dependent) 
parameter characteristic for the impurity system. By 
repeating this calculation for other values of U we ob- 
tain, e.g., 7 (C7 = 4) = 1.25, -y(U = 2) = 1.20, and 
j(U = 0) = 1.0. 

For the BA-LDA the limit on the size of the systems 
which can be treated is much less restraining than for tra- 
ditional methods, since one must diagonalize only a non- 
interacting (Kohn-Sham) Hamiltonian. For small sys- 
tems the accuracy attained is clearly inferior to density- 
matrix rcnormalization group [14] or Quantum Monte 
Carlo [15] methods. This situation parallels the one in 
which DFT finds itself in ab initio calculations: Practical 
applications of ab initio DFT usually do not lead to high 
accuracy. Band structures, for example, are more accu- 
rately calculated using GW [16], and properties of small 
molecules are better obtained from CI [17]. However, 



these are computationally expensive methods that place 
great demands on ones resources and are not easily ap- 
plicable to complex and/or large systems (e.g. molecules 
with more than a few atoms). The power of ab initio 
DFT arises from the relative facility with which it is ap- 
plied to large and complex systems. In the same spirit, 
DFT using generalized LDA's may provide a useful alter- 
native to traditional methods for correlated models with 
a large number of sites. In order to further explore these 
possibilities, we are currently applying the BA-LDA to 
the Mott insulating phase of the 1DHM and to a more 
detailed study of the Friedel oscillations around an im- 
purity. Results will be reported in forthcoming publica- 
tions. 
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FIG. 1. Main figure: Total energy calculated from 
BA-LDA (open squares) and exact diagonalization (full cir- 
cles) versus system size, for a 1DHM with open boundary 
conditions and (7 = 6. For even L we take N = L/2 (so that 
n = N/L — 1/2, corresponding to quarter filling). For odd L 
we take N — (L — l)/2. The error bar at L = 14 represents 
the intrinsic error of the parametrization (2), as estimated 
from the inset. For L — 14 the BA-LDA data agree with the 
exact ones within this error. Inset: Exchange-correlation en- 
ergy of the infinite 1DHM as a function of n for several values 
of U. The full curves represent our parametrization (2) with 
(3), and the symbols represent values obtained from numeri- 
cal solution of the Lieb-Wu BA integral equation for U = 3, 6, 
and 9. For U = and U = oo the parametrization is exact. 

FIG. 2. Main figure: Total energy of the 1DHM with a bi- 
nary site potential and periodic boundary conditions, as de- 
scribed in the main text, for U = 6 and n = 1/2. Open 
squares: BA-LDA, full circles: exact diagonalization. In- 
set: exact xc energy and LDA xc energy, evaluated at ex- 
act and LDA densities, respectively. The absolute errors in 
the LDA xc energy and total energy are different because in 
the LDA the Hartree, external potential and noninteracting 
kinetic components of the total energy are also evaluated at 
the LDA density, not at the exact one. However, the xc en- 
ergy is the dominating source of error in the total energy and 
the relative error in the total energy is considerably smaller 
than that in its xc component. 

FIG. 3. Density distribution for a system with L = 201 
sites, N = 101 electrons, {7 = 6, and a localized impurity with 
Vimp = — 1. Main curve: periodic boundary conditions. Inset: 
open boundary conditions. The density of the same system 
without impurity has been subtracted to display clearly the 
impurity-added contribution. 

FIG. 4. Total energy calculated from BA-LDA for the im- 
purity model of Fig. 3, with open (full symbols) and periodic 
(empty symbols) boundary conditions, plotted as a function 
of the inverse system size. Impurity strength vi = — 1 (upper 
curve), and vi = — 10 (lower curve). 
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